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ABSTRACT 


Three dimensional simulation of oblique hypervelocity impact on orbital debris 
shielding places extreme demands on computer resources. Research to date has shown that 
particle models provide the most accurate and efficient means for computer simulation of 
shield design problems. In order to employ a particle based modeling approach to the wall 
plate impact portion of the shield design problem, it is essential that particle codes be 
augmented to represent strength effects. This report describes augmentation of a 
Lagrangian particle hydrodynamics code developed by the principal investigator, to include 
strength effects, allowing for the entire shield impact problem to be represented using a 
single computer code. 
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CHAPTER Is INTRODUCTION 


This report describes work performed under NASA Grant NAG9-946 aimed at the 
development and validation of a new particle-element method for hypervelocity impact 
simulation. Chapter 2 describes the basic method and illustrates its application with two test 
problems. Chapter 3 describes a series of simulations conducted to validate the method for 
representative Whipple shield orbital debris shielding applications. (Chapters 2 and 3 are 
papers in press in the International Journal of Impact Engineering.) Chapter 4 describes 
additional test problems which simulate published hypervelocity impact experiments. The 
final chapter indicates directions for future work. 


NOTE: References are listed at the end of each individual chapter. 



CHAPTER 2: HYBRID PARTICLE-ELEMENT METHOD 



2-2 


A HYBRID PARTICLE-FINITE ELEMENT METHOD FOR 
HYPERVELOCITY IMPACT SIMULATION 

ERIC P. FAHRENTHOLD and BLAISE A. HORBAN 

Department of Mechanical Engineering, University of Texas, Austin, TX 78712 


Summary— Coupled particle-finite element methods have been suggested as an approach to 
modeling particular impact problems not well suited to simulation with conventional Eulerian, 
Lagrangian, or particle codes. An alternative hybrid particle-finite element technique has been 
developed, in which particles are used to model contact-impact and volumetric deformation 
while finite elements are employed to represent interparticle tension forces and elastic-plastic 
deviatoric deformation. The method has been implemented in a three dimensional code and 
applied to simulate representative hypervelocity impact problems. 


INTRODUCTION 

Particle-based numerical impact models [1] can offer distinct advantages over Eulerian [2] and 
Lagrangian [3] hydrocodes in particular hypervelocity impact applications. An example is the 
design of orbital debris shielding [4, 5], where conventional codes have proven difficult to apply 
[6, 7, 8]. Particle models avoid certain problems with mesh distortion and debris transport which 
have* hindered the effective use of Lagrangian and Eulerian codes in the simulation of three 
dimensional impacts on shielded space structures [9]. However particle methods typically 
incorporate kinematically inexact treatments of material history effects such as plasticity and 
fracture. Recent research efforts have been directed at the formulation of coupled parucle-finite 
element methods for hypervelocity impact simulation [10], as well as a variety of other new 

numerical methods [11, 12]. , . /r ,_ TT . ,, . 

Both particle-in-cell (PIC) methods [13] and smooth particle hydrodynamics (SPH) methods 
[14, 15] employ particles which are actually moving interpolation points. An alternative parucle- 
based modeling methodology developed by Fahrenthold and Koo [16, 17] offers a fully 
Lagrangian, energy-based approach to shock physics simulations. This alternative approach, 
labeled Hamiltonian particle hydrodynamics, avoids the tensile and boundary instabilities 
associated with some smooth particle hydrodynamics formulations [18, 19] and the potentially 
diffusive grid-to-particle mapping schemes characteristic of some particle-in-cell methods. 

In recent work, the particle method of Fahrenthold and Koo has been extended, by coupling the 
aforementioned hydrodynamic particle model to a Lagrangian finite element description of material 
strength in the continuum. The resulting hybrid paiticle-finite element model retains all of the 
features (including general contact-impact effects) of Hamiltonian particle hydrodynamics, while in 
addition accounting for tensile strength, elastic shearing strain, plasticity, and continuum damage 
effects important in the simulation of some hypervelocity impact problems. 

The finite element kinematics used here are similar (not identical) to those employed in existing 
Lagrangian hydrocodes [20], for example DYNA3D [3]. The coordinates of certain nearest 
neighbor particles, identified in the reference configuration, determine finite element nodal 
displacements and hence the local elastic shearing strain, the local plastic strain rate, and 
interparticle tensile forces. Normal and deviatoric continuum damage variables are introduced to 
allow for perforation, fragmentation, and fracture in arbitrary geometry's, while no slideline 
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algorithms [21] are employed. Element failure is based on strictly physical criteria, and is not 
dependent on any smoothing length or mesh-to-particle mapping scheme. Unlike other particle 
based strength models, the present formulation incorporates exact Lagrangian finite strain 
kinematics, which have proven effective in Lagrangian finite element codes in characterizing highly 
strength dependent features of impact dynamics problems. 

A three dimensional, parallel implementation of the hybrid numerical method described here has 
been coded, with specific interest in the application of orbital debris shielding design. The sections 
which follow outline the hybrid formulation, and illustrate the method using two example 
simulations with known experimental results. 


MODELING METHODOLOGY 

The hybrid particle-finite element model described here is formulated using an energy method, 
namely Hamiltonian mechanics [22]. Hence the model formulation procedure differs markedly 
from that used in weighted residual finite element techniques or the finite difference (finite volume) 
method. It differs also from the model formulation procedure used with many particle methods, 
which typically rely in part on energy concepts and in part on the differential balance equations for 
the continuum. As in all discrete Hamiltonian methods, the present model . is formulated by 
assembling kinetic and potential energy expressions for the system, describing the relevant 
constraint equations, and introducing Lagrange multipliers, in order to arrive at the final first order 
state (evolution) equations for the system. Unlike the familiar and purely mechanical Hamiltonian 
models seen often in the literature, the present work employs entropy states to model the thermal 
dynamics of the system. One result is that the energy dissipation expression normally used to 
quantify viscous generalized forces is replaced by a set of nonholonomic constraints on the entropy 
evolution. These and other differences from classical Hamiltonian formulations are discussed in 
more detail by Fahrenthold and Koo [16, 17], who also employ a Hamiltonian methodology to 
formulate a particle model of the SPH type. 

The sections which follow discuss the particle and element kinematics, the stored energy 
expressions for the system, and the nonholonomic constraints, and then apply Hamilton’s 
canonical equations to arrive at a state space description of the impact dynamics problem. The 
particles are used to model kinetic energy effects, contact-impact, and theimomechanical volumetric 
deformation, while the finite elements represent strength effects (interparticle tension, elastic shear, 
and plastic deformation). Contact-impact is modeled using penalty forces similar to those employed 
in DYNA3D. Numerical viscosity is introduced to damp all the elastic modes, and numerical heat 
conduction serves to diffuse shock heating. 


KINEMATICS 

This section provides an overview of the particle and element kinematics. The particles of the 
present model are homogeneously deformed, spherical, Lagrangian control volumes. Hence their 
motion is described completely by a scalar deformation gradient (F) and a center of mass position 

vector (c) for each particle. The material time derivatives (F and c ) of these generalized 
coordinates are gener alize d velocities for the Hamiltonian system. In the reference (undeformed) 
configuration for the modeled system, the particles are arranged in a body-centered cubic packing 

scheme. _ . 

The Lagrangian finite elements used here incorporate nine nodes: they are eight-noded 
hexahedra, with a ninth node located (in the reference configuration) at the element centroid. Eight 
"edge centered" particles define the comers of a hexahedra, while a "body centered" particle locates 
the interior node. In other words, the particle center of mass coordinates are also nodal 
coordinates, for intact (uneroded) elements. Each element is subdivided into six separate five- 
noded subelements, by associating the body centered particle for each element with the six separate 
sets of four particles which define the faces of the hexahedron. The volumes of the subelements are 
used in calculating interparticle tensile forces, while the shear deformation of the hexahedron is 
used to determine the deviatoric strain. The hexahedra are used to describe the following 
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Lagrangian finite strain kinematics for the continuum. The elastic deviatoric strain tensor (E® ) for 
an element is defined by 


E e = E -EP (la) 

where EP is the plastic strain tensor, and the deviatoric total strain is [23] 

E = (1/2) (C - I) ; C =J' 2/3 C; J = [ det(C) ] 1/2 (lb,c,d) 

with J the Jacobian of the hexahedron and C the right Cauchy-Green strain tensor. Since the edge 
centered particles define the hexahedra, the center of mass coordinates for those particles are used 
to calc ulate the strain tensor (C) and element Jacobian (J). As a result the only new generalized 
coordinates (internal state variables) introduced by the elements are the six components of the 

symmetric plastic strain tensor (EP ) for each element. 


KINETIC ENERGY 

The kinetic energy (T) of the system is a function of the particle translation and deformation, 
and takes the form [16, 17] 


T = (1/2) 2 [mW- 1 ptf) 2 U® 2 ]=T(p(i), ) (2a) 

i=l 

where n is the number of particles, m is the particle mass, M is the (constant) particle moment of 
inertia in the reference configuration, p is the particle center of mass momentum, H is the particle 
distributed momentum, and the superscript ,! (i) M denotes the ith particle. The generalized momenta 
are related to the corresponding generalized velocities by 

i(i) ra(iH p<i) ; K0 =m0H H<i) 

(2b,c) 


INTERNAL ENERGY 

The conserved potential in this Hamiltonian model is the system internal energy (U). The 
internal energy is partitioned into three parts: (1) a thermomechanical potential (u, an internal 

energy per unit mass) for each particle which depends on the current particle density (p) and 
entropy per unit mass (s), (2) a deviatoric strain energy which depends on the elastic strain tensor 

(E e ) and a deviatoric damage variable (d, 0 < d < 1), and (3) an interparticle tensile strain energy 
which depends on the particle deformation gradients, the subelement volumes, and a normal 
damage variable (D, 0 < D < 1). Specifically 

U = { 2 m(>) u( i >( p(i>, sCO) } + { 5 (1 - d®) Vj^VO) E e(l) : E e(l) } + 
i=l i=l 

{ 2 2 (1/2) (1 - D®) V®^ j) > 2 } (3) 

i=l j=l 0 
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where ru is the number of elements, V®® denotes the reference volume for element "i", 
c o 

and V e ®i)are the current and reference volume for subelement "j" of element "i", [i® and K® 
o . . 

are the shear modulus and the reference bulk modulus for element "i", is the mass weighted 

average of the particle Jacobians for subelement "j" of element "i' , and the notation < > 
represents the Macauley bracket [ < x > = x for x>0, < x > = 0 for x ^ 0]. In view of the 
aforementioned particle kinematics, 

ft(ij) = ^(i j) (F®)) ; u® = u®( F®, S® ) ; S® = m® s® 

(4a,b,c) 

so that the total particle entropy (S) is a generalized coordinate for the system. Likewise, in view of 
the aforementioned element kinematics, 

E e(i) = E e(i) (c (k) > E P® ) . v e(i,j) = V^J) ( c®) ) (4d,e) 

so that the total internal energy has the functional form 

u = U ( F® , S® , c® , EP (i) , D® . d® ) 

(5a) 

The corresponding generalized conservative forces are 


G® = 


au 

3f® 


; e® = 


au 

as® 


g 


(i) = 


au 

ac® 


(5b,c,d) 


where 0® is the temperature for particle "i", while 


r® 


au ™ au 

aD® ’ 11 ” 3d® 


(5e,f) 

are energy release rates associated with normal and shear damage evolution and 

S® = 2p.® (1 - d®) E e ® =-(l/V£®) 


(5g) 


is a deviatoric stress tensor. 


PLASTICITY MODEL 

The plasticity model used here is a rate independent variation of the isochoric finite strain 
formulation of Fahrenthold and Horban [24]. The flow rule for the incremental plastic strain is 
taken as 

AEP = AX W (6a) 

AX = <x-Y>/{2p. (1-d) [ (1/2) W:W] 1/2 } ; t 2 =(1/2)S:S (6b,c) 


where x is the effective shear stress, Y is the yield stress, and 
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W = CP A' + A' CP 


(6d) 


(6e,f) 


A' = A - (1/3) tr(A) I ; A = S CP + CP S 


CP = I + 2 EP 

In general, the yield stress is taken to be a linear function of J and 9 [25]. Here 


(6g) 


Y = Y(J, 0) = <aJ-l> < 1 - P (0 - 0 O )/(0 m -e 0 )>Y 0 
(7a) 

where 0 O and 0 m are a reference and melting temperature, Y 0 is a reference yield stress, and a 

and P are constants. For use in determining element failure, the effective plastic strain (eP ) is 
tracked by integrating 

eP = [ (1/2) (EP : EP ) ] 1/2 (7b) 

It should be noted that the present modeling methodology allows for the introduction of alternative 
elastic-plastic formulations [26]. 


DAMAGE EVOLUTION EQUATIONS 

The normal and shear damage variables serve to degrade the element moduli in shear and 
tension, and thereby represent the loss of cohesive strength associated with material failure. The 
energy released in damage evolution is a source of irreversible entropy production, so that no 
internal energy is discarded. The simple rate independent damage evolution relations used here are 

AD = Ad = O) ; co = constant (normally 0.1) (8) 

for any time step after an element "fails", due to any one of several effects: (1) the (negative) tensile 
pressure drops below a specified value, (2) the effective shear stress exceeds a specified value, (3) 
the maximum eigenvalue of the deviatoric stress tensor exceeds a specified value, or (4) the 

accumulated plastic strain exceeds a specified value. For to = 0.1 element failure will occur 
gradually, over ten time steps [27]. Once the maximum damage value of 1.0 is reached, the 

element has lost all cohesion and co is again set to zero. 

It should be noted that nothing in the methodology described here precludes the implementation 
of more complex damage evolution models, like those previously implemented in partrcle based 
[28] or Eulerian [29] codes. 


NUMERICAL VISCOSITY 

Considering the preceding discussion of particle and element kinematics, numerical viscosity 
[30] is required to damp both relative motion of the particle mass centers and bulk deformation of 
the individual particles. The damping force on particle "i" due to relative particle motion is taken as 

f(i) = Z vOd) { (c(0 - <$) ) • (c^ - cti)) } (c(0 - c^j) ) / (c® - <$) ) 2 
j=l 


(9a) 
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where is the numerical viscosity 

vW) = C 0 (1/2) ( p® <®A® + p® <!pA® ) A[C (l ’j^ ] (9b) 

;(U) = (hO) +h(i) ) - ic« - c(i) I (9c) 

with c 0 a dimensionless viscosity coefficient, the sound speed for particle "i". A® the cross 

sectional area for particle "i", h® the radius for particle "i", and A[Q a step function ensuring that 
only neighboring particles interact 

A[Q = 1 for C^O ; A[Q = 0 for C<0 (9d) 

Note that the preceding force depends on the normal component of the relative velocity between 
two particles. A similar viscous force which depends on the tangent relative velocity component is 
introduced in each intact element, for relative motion of the edge centered with respect to the body 
centered particle. 

To damp the bulk deformation of individual particles, a viscous pressure is introduced with the 
functional form 

P B 0) =- Cl p® <® h®F® (9e) 

where cj is a dimensionless viscosity coefficient 


NUMERICAL CONDUCTION 

As is standard in impact codes, a numerical heat conduction or artificial viscosity [30] is used in 
the present model. Given the use of entropy states in the Hamiltonian formulation, it takes the form 

n 

jjcon(i) _ (i/Q(i) ) £ R(i,j) (qO) . e(j)) (10a) 

j=l 

where R®j) is the numerical conduction coefficient 

rOJ) = c 2 (1/2) ( p® c® A® + p® c® c® A® ) Aft* 1 ® ] (10b) 

with c® the specific heat for particle "i" and c 2 a dimensionless conduction coefficient 


MECHANICAL AND THERMAL CONSTRAINTS 

The energy balance equation used in conventional continuum codes is replaced here by 
nonholonomic entropy evolution constraints for the particles 

g(0 _ jjirr(i) _ gCon(i) ( lla ) 

where S* rr ® is the entropy production rate due to viscous dissipation, plastic flow, and damage 
evolution 
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s irr(i) = (i/eG) )wG) 


(lib) 


wG) = fG). C® - V® P B G) fG) + I® 6G) + riG) d® + V e G) sG) : EP (11c) 
o o 

with wG) and V® the energy dissipation rate and reference volume for particle "i". As shown by 

Fahrenthold and Koo [16] in the hydrodynamic case, the Lagrange multipliers associated with 
these constraints are known, so that they introduce no new state variables. Similarly, the evolution 
equations for the plastic and damage variables are constraint equations which can be shown to 
introduce no new unknown Lagrange multipliers into the formulation. 

Additional mechanical constraints are must be introduced to represent particle collisions. They 
are 


lc® - c® I - \ (h(i) +h(i) ) > 0 (12a) 

where the constant ^ allows for close packing of the particles at the reference density. The value of 
^ depends on the choice of relative sizes for the edge centered and body centered particles (for 

equal size £ = 0.879). Since the preceding expression is an inequality constraint, it is represented 
in the model formulation process using the nonholonomic expressions 


[ (c® - cCJ) ) / IcO) - c ® I ] • (cG) - c(i)) - \ { h® F® + h® F<J) } = 0 

(12b) 

where h® is the particle radius in the reference configuration. Since the Lagrange multipliers 

associated with the particle collision constraints are numerous, penalty forces A.G® are introduced 
to impose the latter 

JtG j) = kG d) [ £ (hG) +h® ) - lc® - C® i ] A[C c(i d) ] (12c) 


where the step function ensures interaction of overlapping particles only 

£ c GJ) = ^ (h® +h® ) - lc® - c® ] 


(12d) 


and kG J) is a penalty stiffness. For particles connected by springs in series 

kGj) =C 3 / [1.0/(kG) a® 2 / V®) + 1.0/(K(i) A® 2 / V® ] (12e) 


where A^ is the reference cross sectional area for particle "i", and c 3 is a dimensionless penalty 
stiffness, while for parallel springs (as in DYNA3D) 

kGj) = c 3 (1/2) [ kG) a® 2 / v ® + kG) A® 2 / V ® ] 


(120 


HAMILTON'S EQUATIONS 
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The preceding sections have defined the kinematics, stored energy functions, and constraints 
for the physical system. They lead to Hamilton's equations in the form 


(13c) 


p(i) = - g® - f® + L X®j) (c® - c®) / IcCi) - c® I 

j=l 

c® = mOH p® 

if® =-G® -V® P B ® - L Eh® X®# 

j=l ° 


F® =J®-1 H® 


(13a) 

(13b) 


(13d) 


S® = s irr ® - S con(i) (I3e) 

augmented by the evolution equations for the plastic strain and continuum damage variables. 
Integration of these nonlinear, history dependent relations for a chosen equation of state describes 
the thermomechanical dynamics of the impact problem of interest. 

The state space model described here has been implemented in a three dimensional impact code 
[31], developed with a particular interest in orbital debris shielding applications. The state 
equations are integrated using a second order Runge-Kutta method, using time step limits 
described by Fahrenthold and Koo [16]. Currently Mie-Grvmeisen and ideal gas equations of state 
are employed, and linked lists [32] are used to identify neighbor particles. The code incorporates 
compiler directives for parallel execution on Cray systems and SGI workstations. A pre-processor 
is included for model generation, as well as an automated rezoner which deletes particles moving 
outside a user specified control volume. The latter feature is essential for many orbital debris 
shielding simulations. Post processing is performed using commercial graphics software. Testing 
and further code development work is now in progress. The next section presents two example 
simulations, for which there are known experimental results. 


EXAMPLE SIMULATIONS 

The simulations discussed in this section employed material property data from Steinberg [25], 
and the following dimensionless coefficients: 

c Q =c { =0.01 ; c 2 =0.1 ; c 3 =a = p=1.0 ; £ = 0.900 

The material failure stress in shear was set to the maximum yield stress, while the failure pressure 
in tension was set to the spall stress. The effective plastic strain at failure was set to 3.0. 

The first example is an oblique Whipple shield impact simulation, at a velocity of seven 
kilometers per second. The problem parameters are provided in Table 1, and the simulation is 
depicted in Figures 1 and 2. Automated rezoning was used every 1,000 times steps to delete 
particles which move outside the modeled region above the wall plate. For the modeled aluminum 
materials, impact velocity, impact obliquity, shield and wall thickness, and standoff distance, 
experimental ballistic limit curves [33] predict failure of the wall plate for projectiles over 0.45 cm 
in diam eter The modeled projectile of 0.60 cm diameter clearly fails the wail plate, as shown in 
Figure 2, depicting the simulation results at 30.7 microseconds after impact. The required CPU 
time indicated in Table 1 is reasonable, for a three dimensional simulation. Note that the 
requirement to model fragmentation of the projectile, as well as contact-impact of all the projectile 
and shield fragments, presents significant difficulties for conventional Lagrangian codes. On the 
other hand, tracking small debris fragments requires a relatively fine Eulerian mesh. 
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Table 1. Oblique Whipple shield impact simulation 

Projectile diameter (aluminum 6061-T6 sphere) 

= 0.60 cm 

Shield thickness (aluminum 6061-T6) 

= 

0.127 cm 

Wall thickness (aluminum 6061-T6) 


= 0.3175 cm 

Shield-to-wall spacing 

= 

5.0 cm 

Impact velocity 


= 7.0 km/sec 

Impact obliquity 


= 15 degrees 

Number of particles 


= 207,363 

Total simulation time 


= 30.7 microseconds 

Number of time steps 


= 7,000 

Wall clock time (2-CPU SGI Octane) 

— 

60 hours 



Fig. 1. Oblique Whipple shield impact simulation (t = 0.0 microseconds). 



Fig. 2. Oblique Whipple shield impact simulation (t = 30.7 microseconds). 
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Table 2. Long rod impact simulation 

Projectile diameter (DU 0.75% Ti cylinder) 

= 

0.767 cm 

Projectile length 

= 

7.67 cm 

Plate thickness (4340 steel) 


0.64 cm 

Projectile velocity 

= 

0.121 cm/(isec 

Plate velocity 

= 

0.0217 cm/psec 

Impact obliquity 

= 

73.5 degrees 

Number of particles 

= 

66,478 

Total simulation time 

— 

100 (isec 

Number of time steps 

— 

5,934 

Wall clock time (6-CPU SGI Onyx) 

= 

7.50 hours 



Fig. 3. Long rod impact simulation (t = 0.0 microseconds). 



Fig. 4. Long rod impact simulation (Tialf model at t = 100 microseconds). 
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computationally expensive in three dimensions. Hertel [6] reported a requirement of about 50 CPU 
hours on a Cray YMP for a three dimensional Whipple shield simulation using CTH, albeit at a 
standoff distance approximately twice that shown in this example. The latter simulation work 
incorporated significant user intervention, emphasizing the motivation for considering alternative 
methods in the simulation of shielding design problems. 

The second example involves a highly oblique long rod impact on a flat plate, at projectile and 
target velocities of 1.21 and 0.217 kilometers per second. The problem parameters are provided in 
Table 2, and the simulation is depicted in Figures 3 and 4. This problem has been modeled with 
several codes, including CTH [34]. In the present case, pre-processor limitations did not allow for 
the hemispherical nose of the experimental projectile to be modeled, so a flat nose was represented. 
Figure 3 shows the projectile-target configuration at impact, while Figure 4 shows half the physical 
model (cut along the plane of symmetry) at 100 microseconds after impact Experimental data at the 
latter time indicated an eroded rod length of 5.55 cm and a residual velocity of 1.069 kilometers per 
second. The present work yielded a residual rod length of 4.49 cm and a residual velocity of 0.950 
kilometers per second. The corresponding CTH simulation yielded a better estimate of eroded rod 
length (6.00 cm) and residual velocity (0.956 kilometers per second). As in the first example, the 
CPU times reported for a parallel workstation in Table 2 and for a Cray YMP using CTH (4.02 
hours) are not directly comparable. However the CPU time reported here is reasonable, for a 
relatively long (100 microsecond) impact simulation in three dimensions. 


CONCLUSION 

The present paper has outlined the development of a hybrid particle-finite element modeling 
approach to hypervelocity impact simulation. The method appears to have certain advantages oyer 
pure Lagrangian, Eulerian, and particle methods in the application of orbital debris shielding 
design. Although further development and testing of the method is in progress, comparisons of 
simulation results to representative hypervelocity impact problems show reasonable agreement with 
experiment Recent research focused on Arbitrary Lagrangian-Eulerian (ALE), element-free 
Galerkin (EFG), Fluid Implicit Particle (FLIP), and other coupled or hybrid methods suggests that 
such methods will provide opportunities for the expanded use of simulation in the study of 
hypervelocity impact problems. 
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NUMERICAL SIMULATION OF OBLIQUE IMPACT ON ORBITAL 

DEBRIS SHIELDING 

ROBERT J. RABB and ERIC P. FAHRENTHOLD 

Department of Mechanical Engineering, University of Texas, Austin, TX 78712 


Summary — Hybrid particle-finite element methods have been proposed as a modeling 
methodology well suited to the problem of hypervelocity impact simulation. To evaluate the 
use of such numerical methods for orbital debris shielding design, a series of simulations have 
been conducted using a three dimensional hybrid particle-finite element code now under 
development Two sets of oblique impact simulations, one for a single bumper Whipple shield 
and one for a dual bumpier or stuffed Whipple shield, have been compiared to published ballistic 
limit equations, the latter derived from experiment. The results indicate that hybrid particle- 
finite element methods can provide an accurate and computationally tractable approach to the 
simulation of orbital debris shield performance. 


INTRODUCTION 

The design of high performance orbital debris shielding for space structures calls for the 
evaluation of a wide range of new shielding materials and geometry’s. In particular, multi-plate 
geometry’s and composite materials appear to offer significant performance improvements over 
conventional aluminum Whipple shields, using a protection per unit areal weight metric [1]. To 
date most research on the shielding design problem has been pursued experimentally. However 
several factors motivate the development of improved computer aided design tools: (1) promising 
new materials and shield geometry’s have greatly expanded the number of candidate shielding 
designs, (2) the relatively complex interaction of impact velocity, impact obliquity, and material 
failure effects means that a rather large number of experiments are needed to fully characterize the 
three dimensional performance of any single design concept, (3) the desire for faster, cheaper 
spacecraft design places increased emphasis on simulation as opposed to experiment, and (4) the 
range of impact velocities and kinetic energies of interest goes beyond the capabilities of 
conventional light gas guns. 

Although experimental research will continue to occupy a critical role in orbital debris shield 
design, improved computer codes tailored to address the orbital debris impact problem are 
expected to assume a larger role. Experience to date has shown that conventional Eulerian and 
Lagrangian hydrocodes are not ideally suited to this simulation task [2, 3], in particular where 
oblique impact (fully three dimensional) simulations are concerned. Pure particle methods are well 
suited to three dimensional modeling of the debris propagation portion of the shielding design 
problem. However their rather approximate treatment of material strength effects can introduce 
uncertainties in predicting the spacecraft structural response (fracture, spallation, etc.). The latter 
response can be strongly dependent on strength and material history effects, especially where “near 
ballistic limit” simulations are concerned. 

Recognizing the advantages of particle methods in debris propagation calculations and the 
sensitivity of the structural response simulation to material strength effects, some code 
development research has been directed at coupled particle-finite element codes [4] or hybrid 
particle-finite element methods for use in the simulation of hypervelocity impact. The present paper 
describes the application of one hybrid particle-finite element method, described in detail by 
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Fahrenthold and Horban [5], to the specific problem of orbital debris shield design. A series of 
three dimensional simulations involving oblique impact on single and multi-plate aluminum shields 
has been performed for comparison to published experimental ballistic limit curves. The results of 
these simulations demonstrate the utility of particle-finite element methods in the development of a 
new computer codes tailored to the orbital debris shielding design problem. 


NUMERICAL MODELING OF ORBITAL DEBRIS SHIELDING 

Space structures face a current and perhaps growing threat of impact with orbiting space debris. 
For example, inspections of the Space Shuttle orbiters often reveal some impact damage. The 
frequency of these impacts suggests that the population of small orbital debris is increasing, due to 
both spacecraft breakups and collisions among existing orbital debris fragments. For long lived 
space structures, including low earth orbit satellites and the International Space Station (ISS), the 
probability of a serious impact is significant The greatest danger arises from fragments one 
millimeter to one centimeter in length, with a density of 2.8 gm/cc and impact velocities as high as 
13 km/s. Larger fragments exist, but the probability of a collision with debris over one centimeter 
in size is remote [6]. As a result, the focus of much shielding design work, as well as the present 
simulation work, has been on sub-centimeter sized aluminum particle impacts over the 
aforementioned velocity range. 

Most orbital shielding design research has been experimental, with analytical and numerical 
modeling playing a secondary role. There are a number of factors which make numerical 
simulation of oblique impact on orbital debris shielding a very demanding computational problem: 
(1) equation of state and anisotropic strength models are needed for new composite shielding 
materials, (2) better “material failure” models are needed for both composite and metallic shields, 

(3) conventional Lagrangian, Eulerian, and particle codes can be difficult to apply effectively, and 

(4) computer resource requirements are very large. The large memory and CPU time requirements 
are due to: (1) oblique impact (three dimensional) effects, (2) large differences in problem 
characteristic lengths (ratio of shield thickness to standoff distance), (3) a wide range of material 
effects and corresponding time scales (shock thermodynamics govern shield perforation, while 
plasticity and fracture govern structural failure of the waU plate), and (4) fluid-structure interaction 
effects present in pressurized vessel impacts. 

The complexity of the numerical modeling problem can appreciated by considering the simplest 
case of a conventional Whipple shield impact. The simulation problem can be subdivided into three 
phases. The first phase is shield perforation, a “hydrodynamic” event in which material strength 
effects are secondary. The second phase is a debris transport calculation, a non-continuum 
problem which demands very general contact-impact models. The third phase is the wall plate 
impact simulation, where material strength effects are critical and structural failure may occur 
relatively slowly. No single Eulerian, particle-based, or Lagrangian code is ideally suited to 
simulate all three phases of the problem. 

Lagrangian finite element codes normally model perforation by “eroding” elements [7], often 
discarding the internal energy of eroded elements and attaching the eroded element mass to 
neighboring nodes. Their slideline-based contact-impact algorithms can lack the generality needed 
for debris cloud transport calculations, and mesh distortion can mandate frequent rezoning. On the 
positive side, Lagrangian finite element codes provide very accurate material strength models. 

Eulerian finite difference (finite volume) codes model perforation as a hydrodynamic event, 
providing a robust characterization of multi-material impact, at the cost of complex mixed-cell 
thermodynamics and interface tracking algorithms. An inordinately fine mesh may be required to 
model debris transport [8]. Finally, strength modeling is approximate, since in general a different 
collection of material can contribute to the calculation of local material history data at each time 
step. . 

Smooth particle hydrodynamic (SPH) codes [9-12] use “Lagrangian” particles but an “Eulerian” 
internal energy. They provide a very general characterization of contact-impact, although 
conventional SPH density calculations imply a multi-material mixing effect which is usually 
ignored. Particle models are ideal for debris transport calculations, but some problems have been 
experienced with tensile and boundary instabilities [13]. In addition, SPH modeling of strength 
effects is approximate: different material particle sets can contribute to the calculation of local 
material history data at each time step, and numerical “fracture” can occur when particles lose 
contact. 



3-4 


Recognizing the tremendous accomplishments of pure Eulerian (e.g. CTH [14]), Lagrangian 
(e.g. DYNA3D [15]), and particle (e.g. SPHINX [16]) codes, the preceding discussion suggests 
nonetheless that some hybrid formulation may be advantageous for oblique debris shield impact 
simulations. One possibility is some combination of particle and finite element methods - the 
particle capability can simulate shield perforation and debris transport while the finite element 
capability can simulate strength dependent (near ballistic limit) effects during the wall plate impact 
Published work in this area has emphasized adding SPH options to finite element codes [4]. The 
questions which arise with such an approach include: (1) are conventional erosion and slideline 
algorithms sufficiently general? and (2) how will particle-element interactions and element-to- 
particle transition be represented? The answers to these questions will determine in part the range 
of contact-impact problems which can be efficiently modeled using this approach. 

As an alternative approach, another hybrid particle-finite element formulation has been 
developed which combines the Lagrangian particle dynamics method of Fahrenthold and Koo [17, 
18] with a large strain elastic-plastic finite element formulation. This method, described in detail by 
Fahrenthold and Horban [5], has been implemented in a three dimensional computer code (EXOS) 
developed for the simulation of hypervelocity impact on orbital debris shielding. The code can 
simulate, the large strain, elastic-plastic, thermomechanical impact dynamics of solid, fluid, or 
combined solid-fluid structures and systems. The code is particle based, with a finite element 
based strength model, and uses a Mie-Gruneisen equation of state. The kinematics, energy 
functions, and constraints are fully Lagrangian in form. The formulation includes finite strain 
plasticity, scalar damage variables, decoupled volumetric-deviatoric response, and transient 
thermal dynamics. The code is three dimensional and includes a plane of symmetry option. A pre- 
processor simplifies model generation. Although the code is currently undergoing additional 
development, it has been used to perform the simulations discussed in the sections which follow. 

It should be noted here that a number of relatively new numerical modeling techniques including 
Arbitrary Lagrangian-Eulerian (ALE) methods [19], element free Galerkin methods [20], and 
alternative particle methods [21] are under development which may prove effective in addressing 
the orbital debris simulation problem. However to date these methods have not been extensively 
applied to the application of interest here. 




Fig. 1. Example oblique Whipple shield 
impact simulation (t = 0 microseconds). 




Fig. 2. Example oblique Whipple shield 
impact simulation (t = 33 microseconds). 


SIMULATION METHODOLOGY 

As an illustration of the simulations discussed in the sections that follow. Figures 1 and 2 depict 
the results of a typical oblique Whipple shield impact simulation. The parameters of this example 
simulation are provided in Table 1 (material properties were taken from Steinberg [22]). The 
example involves the 7 km/sec impact of a 0.60 cm diameter aluminum sphere on an aluminum 
Whipple shield and wall plate combination, at an impact obliquity of 15 degrees. Figures 1 and 2 
are plots of the simulation results at 0 and 33 microseconds after impact respectively. Consistent 
with the experimentally derived ballistic limit curve of Christiansen [1], the simulation predicts 


perforation of the wall plate. The simulation involved approximately 30,000 particles and required 
approximately 10 CPU hours on a high performance workstation. Automated rezoning was used 
every 500 time steps to minimize CPU time (particles moving outside a user specified region above 
the wall plate were zoned out of the calculation). 


Table 1. Example oblique Whipple shield impact simulation 

Projectile diameter (aluminum sphere) 


0.60 cm 

Shield thickness (aluminum) 

= 

0.127 cm 

Wall thickness (aluminum) 

= 

0.3175 cm 

Shield-to-wall spacing 

= 

5.0 cm 

Impact velocity 

— 

7.0 km/sec 

Impact obliquity 

= 

15 degrees 

Equation of state type 

= 

Mie-Gruneisen 

Spall stress 


0.012 Mbar 

Failure pressure (tension) 

= 

0.012 Mbar 

Yield stress 

= 

0.0029 Mbar 

Plastic failure strain 

— 

3.0 

Shear modulus 

= 

0.271 Mbar 

Density 

= 

2.7 g/cc 

Melt temperature 

= 

1,220 degrees K 

Number of particles 


29,422 

Total simulation time 

= 

33 microseconds 

Number of time steps 

— 

4,500 


AL6061-T6 

Spherical 

Projectile 



* 


0.127 cm 
A1 6061 -T6 
Bumper 


5.08 cm 
Spacing 


AL 6061-T6 

Spherical 

Projectile 


6 


=15° 


0.16 cm 
A1 6061 -T6 
Bumper 


5.715 cm 


11.43 exxu- 




0.32 cm 

A16061-T6 

Bumper 


5.715 cm 


i 0.3175 cm 



0.48 cm 
AL 2219T87 
Rear Wall 


Fig. 3. Geometry for the Whipple shield 
ballistic limit simulations. 


Fig. 4. Geometry for the multi-plate shield 
ballistic limit simulations. 
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In order to conduct a systematic evaluation of the utility of hybrid particle-finite element 
methods in orbital debris shield design, a series of simulations like the one just described were 
performed, in an attempt to match several experimentally derived ballistic limit curves. Such 
ballistic limit curves portray graphically the critical projectile diameter required to fail a specific 
shield and wall plate configuration, as a function of projectile velocity and impact obliquity. Failure 
is defined as perforation or spall of the wall plate. Ballistic limit curves have been constructed by 
testing a wide variety of shield geometry’s and materials, and are used extensively in orbital debris 
shield design. 

The next section discusses the results of a series of Whipple shield (Figure 3) impact 
simulations performed at impact obliquity’s of 15 and 45 degrees, for various projectile diameters 
and impact velocities. The results are compared to known ballistic limit curves derived from 
experiment In the section which then follows, the results of several simulations involving all- 
aluminum stuffed Whipple shields (Figure 4) are compared to corresponding experiments and a 
published ballistic limit curve for a weight equivalent composite-aluminum stuffed Whipple shield. 
The final section notes the conclusions of this study and plans for future work. 


WHIPPLE SHIELD BALLISTIC LIMIT SIMULATIONS 

The Whipple shield is the basis for most simple orbital debris shielding designs. ’Die Whipple 
shield is a thin, sacrificial sheet of material (a bumper) mounted at a fixed standoff distance from 
the protected structure. An incoming projectile hits the bumper and hopefully fragments or melts, 
distributing its impact momentum over a wide area of the rear wall. Experiments on aluminum 
have shown the existence of three distinct projectile response regimes, based on the normal 
component of the projectile impact velocity. The first regime is located below 3 km/s. At low 
velocities, the projectile remains nearly intact after impact with the bumper. Shock pressures are 
low, so the rear wall suffers impact from a deformed but sound projectile. Hence in this regime the 
critical particle diameter decreases as velocity increases. The second (intermediate) velocity regime 
is between 3 and 7 km/s. In this velocity range, the projectile fragments and/or melts on impact, 
reducing its lethality as its velocity increases. The third region consists of normal velocities above 
7 km/s, wherein the debris cloud that impacts the rear wall is a multiphase mixture of projectile and 
bumper components. In this regime the debris cloud becomes more damaging as the impact 
velocity increases. The solid lines in Figures 5 and 6 are experimental ballistic limit curves [1] for 
the Whipple shield configuration of Figure 3, for impact obliquity’s of 15 and 45 degrees 
respectively. 

Numerical impact simulations were performed for the Whipple shield configuration of Figure 3, 
for ten different projectile diameter and impact velocity combinations, at each of the two 
aforementioned obliquity’s, in an attempt to match the experimental ballistic limit curves of Figures 
5 and 6. In particular, simulations were conducted for projectile diameters slightly smaller and 
slightly larger than the critical diameter, at velocities of 3, 5, 7, 9, and 1 1 km/s. In general the 
code was run to a stop time 6-15 times that required for an unimpeded projectile traveling at the 
initial velocity to traverse the space between the shield and rear wall. This allowed sufficient time 
for the debris cloud to impact the rear wall and generate an impulse load. Due to the large number 
of simulations required and the large CPU time requirements of three dimensional simulation, 
relatively coarse models of less than 10,000 particles were initially used. In those cases where the 
coarse models failed to match the ballistic limit curves, finer models (up to 40,000 particles) were 
used in an attempt to improve the simulation results. CPU times for die coarse and fine models 
were approximately 20 and 67 CPU hours respectively, on a Cray J90. Particle count and CPU 
time for each simulation are detailed by Rabb [23]. The simulation results are shown by the open 
and closed triangles plotted in Figures 5 and 6, indicating good agreement, except at die highest 
velocity in the 15 degree obliquity case. Tables 2 and 3 describe the wall plate damage predicted by 
the simulation, including bulging, spallation, and/or perforation. 


MULTI-PLATE SHIELD SIMULATIONS 

One relatively new system designed to improve spacecraft protection from orbital debris is the 
multi-plate or stuffed Whipple shield. The simplest stuffed Whipple shield improves upon the 
basic Whipple by adding an intermediate layer of material between the rear wall and the aluminum 
bumper (Figure 4). This intermediate layer may be aluminum or a composite material, such as 
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0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 

Impact Velocity (km/s) 


Fig. 5. Ballistic limit curve for the 15 degree obliquity Whipple shield simulations. 


Table 2. Predicted damage for the 15 degree obliquity Whipple shield simulations 
(velocity in km/sec, projectile diameter in cm) 


Velocity 

Shield 

Rear Wall Damage Description 

(diameter) 

Failure 


3 (0.20) 

No 

Slight bulge (1 mm deep) 

3 (0.30) 

Yes 

Spallation " 

5 (0.30) 

No 

Slight bulge (1 mm deep) 

5 (0.40) 

Yes 

Spallation; Bulge (2 mm deep) 

7 (0.40) 

No 

Slight bulge (2 mm deep) 

7 (0.50) 

Yes 

Perforation; Hole (17 mm diam) 

9 (0.35) 

No 

Slight bulge (2 mm deep) 

9 (0.45) 

Yes 

Spallation; Bulge (25 mm diam, 9 mm deep) 

11 (0.30) 

Yes 

Perforation; Hole (5 mm diam) 

11 (0.40) 

Yes 

Spallation; Bulge (12 mm diam, 2 mm deep) 
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Fig. 6. Ballistic limit curve for the 45 degree obliquity Whipple shield simulations. 


Table 3. Predicted damage for the 45 degree obliquity Whipple shield simulations 
(velocity in km/sec, projectile diameter in cm) 


Velocity 

(diameter) 

Shield 

Failure 

Rear Wall Damage Description 

3 (0.35) 

No 

Slight bulge (1 mm deep) 

3 (0.45) 

Yes 

Spallation 

5 (0.30) 

No 

Slight bulge (2 mm deep) 

5 (0.40) 

Yes 

Spallation 

7 (0.35) 

No 

Slight bulge (2 mm deep) 

7 (0.45) 

Yes 

Perforation; Bulge (2 mm deep) 

9 (0.40) 

No 

Slight bulge (1 mm deep) 

9 (0.50) 

Yes 

Spallation; Bulge (3 mm deep) 

11 (0.375) 

No 

Slight bulge (1 mm deep) 

11 (0.425) 

Yes 

Spallation; Bulge (18 mm diam, 3 mm deep) 
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combination of Nextel and Kevlar. The intermediate layer can generate additional shocks in the 
projectile and disrupt the moving debris cloud, rendering it less lethal when it reaches the rear wall 
[24, 25]. Although the additional shield does contribute particles to the debris cloud, these 
particles are normally much smaller than the projectile fragments. Experiments have shown that the 
benefits of the intermediate layer outweigh its disadvantages. 

Ballistic limit equations for the stuffed Whipple shield arc not as well developed as those for 
simple Whipple shields. However Christiansen and Kerr [25] do provide ballistic limit equations 
for a specific Nextel- Kevlar-MLI (Multi-Layer Insulation) stuffed Whipple configuration which is 
of potential interest for ISS shielding applications. The solid line in Figure 7 shows that ballistic 
limit curve, for a projectile impact obliquity of 15 degrees. The latter authors also cite several 
experiments on weight equivalent all-aluminum stuffed Whipple shields, demonstrating the 
superior performance of the composite design. Figure 4 shows the weight equivalent all-aluminum 
stuffed Whipple shield used in the cited experiments. 

In the interest of evaluating the particle-finite element modeling methodology used here on a 
stuffed Whipple shield configuration, three simulations were performed for the all-aluminum 
shield configuration of Figure 4, at an impact obliquity of 15 degrees, attempting to duplicate the 
cited experimental results. The results were also compared to the weight equivalent Nextel-Kevlar- 
MLI ballistic limit curve shown in Figure 7. No simulations were performed for Nextel-Kevlar- 
MLI stuffed Whipple designs, due to the lack of equation of state data for the composite materials 
and the fact that the required anisotropic strength models for the composite materials are still under 
development [26]. 

The solid triangles in Figure 7 indicate the projectile diameters and velocities for which 
simulations were conducted. The simulations were run for 61-66 microseconds, approximately 
four times that required for an unimpeded projectile traveling at the initial velocity to traverse the 
space between the first shield and the rear wall. Approximately 35,000 particles were employed, 
with each simulation requiring about 23 CPU hours on a Cray J90. Each of the simulations 
predicted similar results, namely slight bulging and spallation failure of the wall plate. The latter 
results are generally consistent with the experimental data, in that two of the three experiments 
showed wall plate failure. Note from Figure 7 that the tested velocities and projectile diameters 
varied only slightly between the three tests. In addition, the simulations indicate that the 
performance of all-aluminum stuffed Whipple shields is inferior to that of weight-equivalent 
composite designs, given the predicted failures of the wall plate at particle diameters below the 
experimental Nextel-Kevlar-MLI ballistic limit curve shown in Figure 7. It should be noted that 
due to CPU time constraints, the simulations were terminated after the first indication of wall plate 
failure (spallation). The debris cloud of projectile and bumper particles had not in all cases fully 
loaded the rear wall. The rear wall could incur additional damage from subsequent loading of the 
traveling debris. 



Fig. 7. Ballistic limit curve for the 15 degree obliquity multi-plate shield simulations. 
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CONCLUSIONS 

The present paper has described the application of a new particle-finite element modeling 
method in the simulation of hypervelocity impact on orbital debris shielding. Three dimensional 
simulation results were compared to experimentally derived ballistic limit curves for oblique impact 
on Whipple and multi-plate spacecraft shielding, at velocities ranging from 3 to 1 1 km/s. In general 
the simulations show good agreement with the experimental curves, suggesting the further 
development of coupled or hybrid particle-element methods for hypervelocity impact simulation. 

Although the simulation results presented here are encouraging, they represent preliminary 
work involving a new code and a new numerical modeling scheme. Since the present simulations 
were evaluated using a fail/no-fail criterion only, the accurate simulation of failure modes was not 
critically evaluated. In addition, the Mie-Gruneisen equation of state used here is certainly 
approximate at the high end of the impact velocity range considered. Relatively coarse models were 
employed, refined somewhat whenever initial simulation results were not consistent with the 
experimental curves. All of these issues warrant further study. 

Additional code verification and testing is needed, as well as development and implementation 
work on composite equation of state and anisotropic strength models for new shielding designs. 
However the results presented here do support a trend towards the increased use of computer 
simulation in orbital debris shield design. Until experimental capabilities are developed to reach all 
velocities and impact energies of interest, simulations will provide the only estimates of shield 
performance in some hypervelocity impact regimes. 
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EXAMPLE SIMULATIONS 


using the hypervelocity impact code 

EXOS 


Eric P. Fahrenthold, University of Texas at Austin 


The following simulations illustrate application of the hypervelocity impact code 
EXOS (Fahrenthold, 1998), developed at the University of Texas at Austin under the 
support of the Hypervelocity Impact Technology Facility (HIT-F) at NASA Johnson Space 
Center and the State of Texas Advanced Technology Program. 

The first example involves a highly oblique long rod impact on a flat plate, at 
projectile and target velocities of 1.21 and 0.217 kilometers per second. The problem 
parameters are provided in Table 1, and the simulation is depicted in Figures la through Id 
(all figures show half the physical model, cut along the plane of symmetry). This problem 
has been simulated with several codes, including CTH (Hertel, 1992). Figure la shows the 
projectile-target configuration at impact, while Figures lb and lc show the simulation 
results at 50 at 100 microseconds after impact Experimental data at the latter time indicated 
an eroded rod length of 5.55 cm and a residual velocity of 1.069 kilometers per second. 
The present work yielded a residual rod length of 5.4 cm and a residual velocity of 1.0 
kilometers per second. The corresponding CTH simulation (Hertel, 1992) yielded a 
somewhat less accurate estimate of eroded rod length (6.00 cm) and residual velocity 
(0.956 kilometers per second). 

The second example involves oblique impact of an aluminum sphere on a flat plate 
at a velocity of 6.56 kilometers per second. The problem parameters are provided in Table 
2, and the simulation is depicted in Figures 2a through 2d. Note that unlike Eulerian and 
SPH (Smooth Particle Hydrodynamics) models, the present true Lagrangian formulation 
makes an unambiguous prediction regarding the hole size in the plate (Figure 2c). An areal 
density plot produced from the simulation (Figure 2d) shows good agreement with the 
radiograph of the corresponding experiment, the latter reported by Piekutowski (1996). 

The third example involves a pair of oblique Whipple shield impact simulations, at a 
velocity of 7.0 kilometers per second. The problem parameters are provided in Table 3, and 
the simulations are depicted in Figures 3a through 3d. For the modeled aluminum 
materials, impact velocity, impact obliquity, shield and wall thickness, and standoff 
distance, experimental ballistic limit curves (Christiansen, 1993) predict failure of the wall 
plate for projectiles over about 0.50 cm in diameter. The modeled projectile of 0.40 cm 
diameter deforms but does not perforate the wall plate, as shown in Figure 3b, depicting 
the simulation results at 54 microseconds after impact. On the other hand, the modeled 
projectile of 0.60 cm diameter clearly fails the wall plate, as shown in Figure 3d, again at 
54 microseconds after impact. Note that the requirement to model fragmentation of the 
projectile, as well as contact-impact of all the projectile and shield fragments, presents 
significant difficulties for conventional Lagrangian codes. On the other hand, tracking small 
debris fragments requires a very fine Eulerian mesh, computationally very expensive in 
three dimensions (Hertel, 1993). The latter simulation work also incorporated significant 
user intervention, emphasizing the motivation for considering alternative methods in the 
simulation of shielding design problems. With regards to SPH methods, poor strength 
modeling generally makes accurate characterization of the wall plate ballistic limit extremely 
difficult (Faraud et al., 1998). It appears that no EFG (Element Free Galerkin) based work 
has attempted to simulate an impact and penetration problem as complex as this example 
(Belytschko et al., 1996). 




Figure lb. Oblique long rod impact, color 
on temperature (t = 50 microseconds). 


Figure Id. Oblique long rod impact, color 
on velocity (t = 100 microseconds). 


Projectile diameter (DU 0.75% Ti cylinder) 

Projectile length 

Plate thickness (4340 steel) 

Projectile velocity 
Plate velocity 
Impact obliquity 
Equation of state 
Number of particles 
Total simulation time 


0.767 cm 
7.67 cm 
0.64 cm 
1.21 km/sec 
0.217 km/sec 
73.5 degrees 
Mie-Gruneisen 
66,478 

100 microseconds 





Figure 2a. Oblique sphere impact, color 
on material (t = 0.0 microseconds). 


Figure 2c. Plate perforation, color on 
plastic strain (t = 6.6 microseconds). 



Figure 2b. Oblique sphere impact, color 
on material (t = 6.6 microseconds). 
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Figure 2d. Oblique sphere impact, 
greyscale on areal density 
(t = 6.6 microseconds). 


Table 2: Oblique plate impact simulation 

Projectile diameter (aluminum sphere) 

Plate thickness (aluminum) 

Impact velocity 
Impact obliquity 
Equation of state 
Number of particles 
Total simulation time 


0.953 cm 
0.1143 cm 
6.56 km/sec 
45 degrees 
SESAME 37 18 
33,146 

6.6 microseconds 








Figure 3a. Oblique Whipple shield Figure 3c. Oblique Whipple shield 

impact, projectile below ballistic limit size impact, projectile above ballistic limit size 

(t = 0 microseconds). (t = 0 microseconds). 







Figure 3b. Oblique Whipple shield 
impact, projectile below ballistic limit size 
(t = 54 microseconds). 


Figure 3d. Oblique Whipple shield 
impact, projectile above ballistic limit size 
(t = 54 microseconds). 


Table 3: Oblique Whipple shield impact simulations 


Projectile diameters (aluminum spheres) 
Shield thickness (aluminum) 

Wall thickness (aluminum) 

Shield to wall spacing 
Impact velocity 
Impact obliquity 
Equation of state 
Number of particles 
Total simulation time 


0.40 & 0.60 cm 
0.127 cm 
0.3175 cm 

5.0 cm 

7.0 km/sec 
15 degrees 
SESAME 3718 

109.000 

54 microseconds 
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The fourth example represents a simulation of Southwest Research Institute test 
number 7139-6, which involved a multi-plate shield. The problem parameters are provided 
in Table 4, and the simulation is depicted in Figures 4a through 4c. The figures show the 
simulation results at 74.0 microseconds after impact 



Projectile mass (hollow cylinder) 
Impact velocity 
Impact obliquity 
Shield thickness (aluminum) 
Wall #1 thickness (aluminum) 
Wall #2 thickness (aluminum) 
Shield- to- wall #1 spacing 
Wall #l-to-wall #2 spacing 
Equation of state type 
Failure strain 
Number of particles 
Total simulation time 


I. 06g 

II. 18 km/sec 
65 degrees 
0.16002 cm 
0.3175 cm 
0.2032 cm 
6.096 cm 
2.53 cm 
SESAME 3718 
1.00 
216,106 

74.0 microseconds 


The simulation shows no wall plate failure, as observed in the experiment Hole 
sizes may be compared as follows (all dimensions in centimeters): 


Shield #1: 
Shield #2: 
Wall plate: 



4.5 x 3.0 (approximate) 
3.0 x 3.0 (approximate) 
bulge, no failure 


The preceding results could perhaps be improved by reducing the strength of the 
inhibited shaped charge projectile, which was modeled here as a hollow sphere of intact 
aluminum. 




Figure 4a. 
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Figure 4b. 
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Figure 4c 
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CHAPTER 5: CONCLUSION 


This report has described work performed under NASA Grant NAG9-946, 
including the development and validation of a new particle-element method for 
hypervelocity impact simulation. Future work is aimed at extension of the method to 
include more complex materials and more complex shielding designs. 



